Volatility-dependent damping of evaporation- driven Benard-Marangoni instability 
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The interface between a pure liquid and its vapor is usually close to saturation temperature, hence 
strongly hindering any thermocapillary flow. In contrast, when the gas phase contains an inert gas 
such as air, surface-tension-driven convection is easily observed. We here reconcile these two facts by 
studying the corresponding crossover experimentally, as a function of a new dimensionless number 
quantifying the degree of damping of interfacial temperature fluctuations. Critical conditions are in 
convincing agreement with a simple nonlocal one-sided model, in quite a range of evaporation rates. 
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When a layer of volatile liquid is exposed to a dry non- 
condensable (or inert) gas, evaporation occurs and the 
heat needed for the phase change induces a cooling of the 
liquid/gas interface. At large enough thickness and/or 
evaporation rate (i.e. large enough Marangoni number 
Ma), surface-tension-driven convection is observed |l|, 
often in the form of hexagonal patterns with their typ- 
ical defects Q. However, the analogy with the classical 
Benard-Marangoni (BM) convection is far from straight- 
forward, given a number of complicating issues partly 
discussed hereafter (see also Q)- Among these, satura- 
tion of the gas phase should play an essential role, given 
that surface temperature variations are a priori excluded 
in a one-component liquid-vapor system (indeed, except 
for extremely small liquid depths, the interface should 
be close to local equilibrium [2[). Another extreme is the 
non- volatile case, where experiments 0] have confirmed 
Pearson's theoretical value of Ma c c± 80 for the criti- 
cal Marangoni number [5[. Hence, Ma c should somehow 
depend on the liquid volatility, a convenient measure of 
which is the mole fraction uo^ = Ps&t/Pg, i- e - the ratio of 
the saturation pressure p sat at ambient temperature to 
the total gas pressure p g . Studying the crossover between 
the nonvolatile regime (ujy: <C 1) and the pure vapor case 
(c<j£ — >- 1) is the main goal of the present Letter. 

In addition to its fundamental interest, pattern forma- 
tion in drying liquid layers or thin films presently lacks 
quantitative understanding, which is however crucial for 
many techniques such as polymer coating [6[ or nanopar- 
ticle deposition Q. An essential step towards the anal- 
ysis of the wide variety of obtained patterns is a correct 
evaluation of the supercriticality e = (Ma — Ma c )/Ma c , 
which clearly requires an accurate determination of Ma c . 
Another challenging objective is to understand the rich 
nonlinear dynamics of evaporation-driven BM patterns 
in order to control and even optimize deposition struc- 
tures, which requires developing simple models of highly 
supercritical convection. This is the second aim of this 
Letter, which presently focuses on the case of a pure liq- 
uid, where convection is driven by thermal effects only. 



In order to understand the mechanism responsible for 
the damping of interfacial temperature fluctuations, as a 
function of the volatility, let us consider a flat interface 
(at z = 0, with z pointing to the gas) at temperature Tu, 
where evaporation occurs at a mass flux density J (in 
kg/m 2 s). Hereafter, a subscript E will denote a quantity 
evaluated at the interface, the gas mixture is taken to be 
perfect, its total pressure p g is supposed to be constant 
and uniform (small dynamic viscosity), and the inert gas 
(say, air) is not absorbed into the liquid. This implies 
(see e.g. Q) 
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where D is the vapor- air diffusion coefficient, M v is the 
molar mass of the vapor, p v is the partial pressure of va- 
por in the gas phase, uo = p v /p g its mole fraction, and R 
is the ideal gas constant. In addition, the energy balance 
at the interface reads 
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where T\ and T g are respectively the liquid and gas tem- 
peratures, C is the latent heat of vaporization, while A/ 
and are respectively the liquid and gas thermal con- 
ductivities (with A^ <C A^ in general). 

We now consider fluctuations (denoted by tilded quan- 
tities) around a particular steady (or quasi-steady) state 
distinguished by a superscript 0. The determination of 
this particular state needs not be detailed for the mo- 
ment, and in principle, the following reasoning applies to 
both evaporation (J° > 0) and condensation (J° < 0). 
The interface temperature is written = T§ + , and 
assuming local chemical equilibrium at the interface, the 
corresponding fluctuation of vapor partial pressure there 
reads p v ^ = Psat(^s) ^£0 where p s &t(T) is the coexistence 
(i.e. Clausius-Clapeyron) curve and a prime denotes its 
derivative. As fluctuations satisfy V 2 p v = in the limit 
of a small Peclet number (defined on a typical length scale 
of the fluctuations, assumed to be much smaller than 
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the typical size H of the gas phase) and in the quasi- 
static hypothesis, we have p VCL = Pg at (T^) T^ q e - ^' 2 , 
where a subscript q indicates the Fourier component with 
wavevector q (in the horizontal plane). Similarly, one 
also has V 2 T g = 0, because the Lewis number Le = D / ' n g 
is 0(1) in the gas (n g is the gas thermal diffusivity ) . 
Hence, assuming T g = T\ (= Te) at z = 0, we have 
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Finally, linearizing Eq. (pQ), we can calculate (the 
Fourier transform of) the phase change rate fluctuation 
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where fluctuations of the denominator have been ne- 
glected (this is rigorously valid for |q| _1 <C H, as will 
be shown elsewhere). Then, Fourier transforming the in- 
ter facial energy balance (|2j) and grouping terms, we get 
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Equation ((4]) has the form of a Newton's cooling law for 
liquid temperature fluctuations, with a heat transfer coef- 
ficient depending upon their wavenumber |q| (hence, the 
physical space expression of Eq. (j4]) involves a nonlocal 
convolution term) . The positive dimensionless number a 
turns out to be an effective gas-to-liquid ratio of thermal 
conductivities, accounting for the effect of phase change 
through its second term. In particular, the latter contri- 
bution is seen to diverge for uj^ —> 1, i.e. in the limit of a 
pure vapor phase, for which Eq. (j4]) implies T^ q = 0, i.e. 
the interface temperature does not fluctuate and remains 
equal to the saturation (boiling) temperature. 

Now, applying Eqs (|4j) and (|5]) to the modeling of 
evaporation-driven BM convection in a liquid layer of 
height e much thinner than the gas phase thickness i/, 
it turns out that Pearson's theory [5[ can be straightfor- 
wardly applied, using an effective Biot number Bi = ak 
where k = |q| e is the dimensionless wavenumber. The 
neutral stability threshold is then given by 

_ 16fc(fccosh[fc] + ak sinh[fc])(2fc - sinh[2fc]) 
ak ~ 4k 3 cosh[/c] + 3 sinh[/c] - sinh[3/c] ( ^ 

The critical Marangoni number Ma c (a) and the critical 
wavenumber k c (a) are then found by minimizing Ma^ 
with respect to fc, and will now be compared to exper- 
iments. Note finally that the theory just described ap- 
pears as a particular case of a more general formulation 
(not limited to e <C H) described in [3|, and also based 
on a one-sided reduction of the evaporation-driven BM 
problem by adiabatic slaving of gas phase fluctuations. 



In order to test these predictions, accurate experiments 
were performed by evaporating thin liquid layers into am- 
bient air (T w 24°C) at rest, until the liquid completely 
disappears. As explained hereafter, we mostly focus on 
the moment at which convective patterns disappear in fa- 
vor of a uniform evaporative state. Volatile liquids used 
are Hydrofluoroethers, HFE-7000, 7100, 7200 and 7300 
from the company 3M, which have similar physical prop- 
erties apart for their saturation pressure p sat (factor of 
about 2 between two successive HFEs). HFE-7000 is the 
most volatile with p sat (24°C) = 0.61 bar and HFE-7300 
is the less volatile with p sat (24°C) = 0.06 bar. Other 
thermodynamic and transport properties used hereafter 
are found in 3M data sheets (available on 3M web site). 

Each experimental run is started by pouring a certain 
amount of HFE in a cylindrical container to form an ap- 
proximately 1 mm thick liquid layer. The container is 
made of a PVC cylinder glued by silicone to a 10 mm 
thick aluminum plate. The height of the cylinder is 1 
cm, its diameter is 63.5 mm and its thickness is 6 mm. 
In addition to the effect of volatility (dependent on the 
HFE used), we also vary the evaporation rate indepen- 
dently by changing the "transfer distance" H in the gas. 
This is accomplished by topping another PVC cylinder 
(of the same diameter) on the one glued to the plate, 
wrapping them with a scotch tape in order to avoid any 
vapor leak. Using additional cylinders of various heights 
allows to set H to 1 cm, 2 cm, 3 cm, 4 cm and 5 cm. 

In these conditions, the evaporation process is limited 
by diffusion of vapor into air and the evaporation rate 
E = J°S (where S is the container cross section) remains 
quasi-constant until the layer is too thin and dewetting 
begins. The liquid film thickness, e, is measured by 
weighting and is deduced from the measured total mass, 
mtot, taking into account the mass of the liquid meniscus 
against the internal cylinder wall, m men , and the mass 
of the vapor contained in the gas phase above the liquid, 
rn v , such that e = (m tot — rri rnen — m v )/piS, where p\ is 
the liquid density. Both these contributions cannot be ne- 
glected because the critical layer thicknesses are generally 
small. More precisely, m men is theoretically estimated as- 
suming that the meniscus is in its hydrostatic equilibrium 
state and that the liquid is perfectly wetting. In turn, m v 
has been measured experimentally for each H and each 
HFE (molecular weight between 200 and 350 g/mol) us- 
ing a suspended thin circular dish filled with liquid and 
placed very close to the container bottom but without 
touching the container wall, hence "simulating" the pres- 
ence of a liquid layer. In the worse case (most volatile 
liquid HFE-7000 and highest container H = 5 cm), we 
find m v ^ 0.2 m to t and in the best one (HFE-7300 and 
H = 1 cm), we get m v ^ 0.005 m to t- The relative un- 
certainty on the liquid layer thickness measurement is 
estimated to be lower than 2.5%. The evaporation rate, 
E, is simply computed from the time derivative of the 
total mass (E = —dm to t/dt) using a linear fit. 
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As convection in the pure liquid is necessarily associ- 
ated with temperature variations, we use a Focal Plane 
Array IR camera-type (Thermosensorik, InSb 640 SM) 
facing the liquid/gas interface, to follow the time evolu- 
tion of the whole cellular pattern. IR images and mtot(t) 
are recorded at a frequency of 1 Hz during the drying 
of the liquid layer. Typically, the observed sequence is 
similar to the one obtained in [1], i.e. convection ap- 
pears right after filling and the pattern is strongly time- 
dependent, evolving into more stable hexagonal-like ar- 
rangements when e decreases, until the convective state 
turns into a "conductive" one. The convection cells do 
not disappear altogether, however. At a certain mo- 
ment, a straight front separating convective and conduc- 
tive states starts to propagate along a horizontal direc- 
tion (at a constant velocity) until convection completely 
disappears. Performing specific experiments in which 
the container tilt angle was intentionally slightly varied 
showed that this front is merely due to a non-absolute 
horizontality of the layer (the front velocity decreases 
when increasing the tilt angle and vice versa). We have 
chosen to define the critical liquid layer thickness, e c , by 
e c = (ei + e2)/2 where e\ is the measured thickness when 
the front starts to propagate and 62 the measured thick- 
ness when convection has totally disappeared. For the 
small tilt angles tested, e c is found to be independent of 
these small deviations w.r.t. absolute horizontality. 

From the measurement of e c and E, we then estimate 
the critical temperature difference across the liquid layer, 
AT C , using the thermal balance ((2]) with a linear tem- 
perature profile in the liquid and neglecting heat com- 
ing from the gas phase, such that AT C = ECe c / XiS. 
Then, the critical Marangoni number Ma c is calculated 
as Ma c = —^T^T c e c /rjiKi, where is the surface ten- 
sion variation with temperature, rji is the liquid dynamic 
viscosity and ki is its thermal diffusivity. The value of 
7t has been measured for each HFE by the pendant drop 
method using the tensiometer Kriiss DSA100 with its 
thermostated chamber, and a thermocouple placed at the 
syringe tip end to measure the drop temperature accu- 
rately. Surface tension has been measured in the range 
15-30°C, taking special care in order to maintain the drop 
in a saturated environment (procedure validated by mea- 
suring 7x of ethanol). 

According to our model, the value found for Ma c 
should only depend on the liquid used (through the value 
of a characterizing the damping of thermal fluctuations 
at the interface) and not on its evaporation rate E. More 
precisely, as AT C ~ E e c for a given liquid, Ma c ~ E e 2 c 
should be a constant, leading to the scaling e c ~ E~ x l 2 . 
Apart for some small variations studied hereafter, the size 
of convection cells at threshold should be roughly propor- 
tional to the depth e c . Hence, the critical wavenumber 
q c ~ e~ l ~ E x l 2 . Both these scalings indeed match ex- 
perimental measurements, as shown in Fig. HJ 

Now, in order to fully validate our Pearson-like theory, 
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FIG. 1. Measured critical liquid thickness e c (left) and critical 
wavenumber q c (right) as a function of the evaporation rate 
E, for different liquids (triangles : HFE- 7000, squares : HFE- 
7100, circles : HFE-7200, crosses : HFE- 7300). Straight lines 
indicate the theoretical scaling laws (see text). 

a is directly computed from Eq. ([5]) using our own mea- 
sured values of D for each HFE, obtained by the Stefan's 
tube method The obtained values of Ma c for all the 
HFEs and for all the container heights H are plotted as a 
function of a in Fig. O and compared to the theoretical 
law Ma c (a), independent of E as already mentioned. 
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FIG. 2. Measured critical Marangoni numbers (symbols) as 
a function of a, for the four different liquids and the five 
container heights. Inset: corresponding measured critical 
wavenumbers. Theoretical laws are shown as plain curves. 

The inset of Fig. [2]also shows the dimensionless critical 
wavenumber of the pattern, k c = q c e c ^ clearly indepen- 
dent of E as well. Note that the critical wavenumber q c 
(as also plotted in Fig. [1]) and the wavenumber q (during 
the convection regime) are measured as the mean po- 
sition of the fundamental peak in azimuthally averaged 
FFT spectra of IR images (see Fig. [3]). Error bars on 
k c in Fig. [2] correspond to the width at middle height 
of the fundamental peak, hence indicative of the level of 
disorder in the pattern (higher at high volatility). 

Figure [2] demonstrates the strong stabilizing effect of 
the liquid volatility, as well as a quite satisfactory agree- 
ment with our simple one-sided theory (given typical un- 
certainties remaining on some fluid properties and the 



absence of fitting parameters). Note that this actually 
validates a number of assumptions made in such type of 
models (see also [3]), such as small gas viscous stresses, 
low Peclet numbers in both phases, quasi-steadiness of 
the approach despite the continuously decreasing liquid 
depth, undeformable interface, absence of temperature 
discontinuity, ... In addition, we emphasize that the sim- 
plest form of the theory presented here relies on the ad- 
ditional assumption of a large gas thickness H compared 
to the liquid depth e. As the length scale of connec- 
tive fluctuations typically scales with e, their penetration 
depth in the gas is of the same order, which in fact al- 
lows neglecting the effects of gas density variations and 
diffusion-induced convection (even though both these ef- 
fects do affect the homogeneous evaporation state, hence 
J°, when uj^ increases). This will be detailed elsewhere. 

To conclude, let us briefly explore nonlinear regimes of 
evaporative BM convection. Figure [3] shows that cellular 
patterns become more regular when the supercriticality 
e = (Ma — Ma c )/Ma c decreases (i.e. when time goes on), 
and that at the same value of e, patterns are more disor- 
dered for more volatile liquids. This is also confirmed by 
the corresponding Fourier spectra, from which averaged 
wavenumbers k = qe depicted in Fig. |4] are extracted. 
We first note that the shape of fc(e) curves (including the 
non-monotonic behavior at low a) is in nice qualitative 
agreement with direct simulations of [9[, which however 
rely on a constant Biot number instead of Eq. (jlj. 




FIG. 3. Typical evaporation- driven BM patterns (H = 3 cm): 
(a) HFE-7000, e = 0.2; (b) HFE-7000, e = 2; (c) HFE-7300, 
e = 0.2; (d) HFE-7300, e = 2. White bars are 4 mm long. 
Top-right insets : contourlines of power (Fourier) spectrum. 

Interestingly, Fig. 2] also shows that the measured 
wavenumber evolutions are rather independent of the 
container height (hence on the evaporation rate), while 
they do depend on the liquid used. This clearly has to do 
with the fact that the timescale for liquid depth variation, 
r e ~ e/|e|, is always much larger than the thermal time 
scale Tth ~ £ 2 /ki (quasi-steady evolution). However, the 
fact that r e turns out to be much smaller than the lat- 



eral diffusion time tl ~ L 2 / ' k\ (L is the container size) 
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FIG. 4. Measured dimensionless wavenumber k versus super- 
criticality e = (Ma — Ma c )/Ma c . For each liquid, evolutions 
corresponding to all five container heights are represented. Pr 
is the Prandtl number of the liquids. 

points to a rather fast mechanism of wavelength selection, 
which might be due, at least far from threshold where the 
pattern is large-scale (k <C fc c ), to the anomalous dissi- 
pation mechanism described by Eq. (j4]). This remains 
to be studied however, along with the quite unexplored 
scenarii of transition to "interfacial turbulence" (see also 
Q), for which the one-sided model we propose should 
be appropriate in a quantitative sense. Note finally that 
although validated here using a Benard set-up, the phase- 
change-induced homogenization mechanism described by 
Eq. (|4]) is expected to be generic for other geometries 
as well (e.g. drops, bubbles, ...), at least for sufficiently 
short-scale interfacial temperature fluctuations. 
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